function [Xi,Ti,Ui,cell_type] = interpolate(grid,base,k,U,varname)
% [Xi,Ti,Ui,cell_type] = interpolate(grid,base,k,U,varname)
%
% Interpolate on a set of Lagrangian nodes. Input and output are
% assumed to be in element base ordering: this allows to treat
% discontinuous functions.
%  varname: can be {'lambda','phi','omega'}
% The results are as follows:
%  Xi: nodal coordinates
%  Ti: connectivity of the reconstructed field
%  Ui: interpolated values
%  cell_type: paraview code for the elements used in the reconstructed
%             field
%
% Note: this function relies on a file localLagmeshes.octxt containing
% the local grids, with local connectivity and cell_type indication.
% The local connectivity can be done by any number of P1 or P2
% elements, which are the two type of elements used by paraview.
%
% Note: this function works for a generic dimension grid.d (as far as
% the corresponding data are provided by localLagmeshes.octxt). No
% fictitious components are introduced (for instance, a third zero
% component for 2D data).

 % 1) Get the nodes
 locmeshes = load('localLagmeshes.octxt');
 if(k>size(locmeshes.meshes,2))
   error(['Degree ',num2str(k),' not supported.']);
 endif
 switch varname

  case "lambda"
   locmesh = locmeshes.meshes(base.me.d-1,k);

  case {"phi","omega"}
   locmesh = locmeshes.meshes(base.me.d,k);

  otherwise
   error (["Unknown varname value ",varname]);
 endswitch

 p = locmesh.p;
 t = locmesh.t;
 cell_type = locmesh.cell_type;

 % 2) Evaluate the basis on the local mesh
 np_tot = size(p,2); % number of points per old cell
 np     = size(t,1); % number of points per new cell
 nt     = size(t,2); % number of new cells per old cell
 switch varname

  case "lambda"
   PHI = zeros(base.nk,np_tot);
   for i=1:base.nk
     PHI(i,:) = ev_pol(base.e_s{i},p);
   end

  case "phi"
   PHI = zeros(base.pk,np_tot);
   for i=1:base.pk
     PHI(i,:) = ev_pol(base.p_s{i},p);
   end

  case "omega"
   PHI = zeros(base.me.d,base.mk,np_tot);
   for id=1:base.me.d
     for i=1:base.mk
       PHI(id,i,:) = permute(ev_pol(base.o_s{id,i},p),[3,1,2]);
     end
   end

 endswitch

 % 3) Build nodes, connectivity and nodal values
 idx = ones(1,np);
 
 switch varname

  case "lambda"
   Xi = zeros(grid.d,np,nt*grid.ns);
   Ti = zeros(       np,nt*grid.ns);
   Ui = zeros(np,nt*grid.ns);
   for is=1:grid.ns
     for it=1:nt
       
       i = (is-1)*nt + it;
       Xi(:,:,i) = grid.s(is).b * p(:,t(:,it)) + grid.s(is).x0(:,idx);

       Ti(:,i) = (is-1)*(nt*np) + (it-1)*np + [1:np];

       Ui(:,i) = PHI(:,t(:,it))' * U(:,is);

     end
   end

  case {"phi","omega"}
   Xi = zeros(grid.d,np,nt*grid.ne);
   Ti = zeros(       np,nt*grid.ne);
   switch varname
    case "phi"
     Ui = zeros(np,nt*grid.ne);
    case "omega"
     Ui = zeros(base.me.d,np,nt*grid.ne);
   endswitch
   for ie=1:grid.ne
     piola = grid.e(ie).b/grid.e(ie).det_b;
     for it=1:nt

       i = (ie-1)*nt + it;
       Xi(:,:,i) = grid.e(ie).b * p(:,t(:,it)) + grid.e(ie).x0(:,idx);

       Ti(:,i) = (ie-1)*(nt*np) + (it-1)*np + [1:np];

       switch varname
        case "phi"
         Ui(:,i) = PHI(:,t(:,it))' * U(:,ie);
        case "omega"
         % include the Piola transformation
         for id=1:base.me.d
           Ui(id,:,i) = (permute(PHI(id,:,t(:,it)),[2,3,1])' * U(:,ie))';
         end
         Ui(:,:,i) = piola*Ui(:,:,i);
       endswitch
     end
   end

 endswitch

return

